%Gitee Repo

clc
clear

n = 1000;
dx = 0.01;
dt = 0.001;
c=1;
T=1;
w=2*pi/T;

sigma = (c*dt/dx)^2;

u0 = zeros(n,1);
u1 = zeros(n,1);
u2 = zeros(n,1);
_x = (linspace(0,dx*(n-1),n))';

figure()
k=0;
imgindex = 0;

for t = 0:dt:10
 %  for i = 2:n-1 %Numerical PDE Iteration
 %    u2(i) = 2*u1(i) - u0(i) + sigma*(u1(i-1)-2*u1(i) + u1(i+1));
 %  endfor

  % 向量化
  u2(2:n-1) = 2*u1(2:n-1) - u0(2:n-1) + sigma*(u1(1:n-2)-2*u1(2:n-1) + u1(3:n));

  if t < T %Boundary Conditions
    u2(1) = sin(w*t);
  else
    u2(1) = 0;
  end
  u2(n)= 0;

  u0 = u1; %Update Time
  u1 = u2;

  if mod(k,50) == 0
    clf
    hold on
    axis equal
    axis([0 10 -2 2])

    plot(_x,u2)
    pause(0.01)
    k;
    imgindex++;
  end

  k++;
end
